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Abstract 

Functional  imaging  of  the  resting  brain  consistently  reveals  broad  motifs  of  correlated  blood 
oxygen  level  dependent  (BOLD)  activity  that  engage  cerebral  regions  from  distinct  functional 
systems.  Yet,  the  neurophysiological  processes  underlying  these  organized,  large-scale 
fluctuations  remain  to  be  uncovered.  Using  magnetoencephalography  (MEG)  imaging  during  rest 
in  12  healthy  subjects  we  analyse  the  resting  state  networks  and  their  underlying  neurophysiology. 
We  first  demonstrate  non-invasively  that  cortical  occurrences  of  high-frequency  oscillatory 
activity  are  conditioned  to  the  phase  of  slower  spontaneous  fluctuations  in  neural  ensembles.  We 
further  show  that  resting-state  networks  emerge  from  synchronized  phase-amplitude  coupling 
across  the  brain.  Overall,  these  findings  suggest  a  unified  principle  of  local-to-global  neural 
signaling  for  long-range  brain  communication. 
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1  Introduction 

Over  the  past  15  years,  there  has  been  considerable  interest  in  studying  the  resting  activity  of 
the  human  brain,  particularly  with  functional  magnetic  resonance  imaging  (fMRI)  (Raichle, 
2011).  One  remarkable  property  of  spontaneous  cerebral  activity  is  that  it  consistently 
segregates  in  resting-state  networks  (RSN),  which  coincide  anatomically  with  the  major 
functional  systems  of  the  brain  (Biswal  et  al.,  2010).  Consequently,  studies  suggest  that  RSN 
provide  insight  into  the  large-scale  mechanisms  of  healthy  and  impaired  neural 
communication  (Buckner  et  al.,  2005;  Tomasi  and  Volkow,  2012). 

Yet,  our  current  understanding  of  the  neurophysiological  mechanisms  underlying  these 
large-scale  fluctuations  remains  elusive  (Leopold  and  Maier,  2012;  Smith,  2012). 
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Simultaneous  BOLD  and  electrophysiology  recordings  in  the  visual  cortex  of  monkeys 
report  that  gamma  activity  of  the  local  field  potentials  (LFP)  is  the  main  electrophysiological 
basis  of  the  BOLD  signal  (Logothetis  et  al.,  2001;  Shmuel  and  Leopold,  2008).  These  local 
results  cannot  resolve  how  the  large-scale  connectivity  across  the  whole  brain  of  the  resting 
state  is  generated.  Human  magnetoencephalography  (MEG)  of  the  whole  brain  emphasized 
the  contribution  of  alpha  and  beta  oscillatory  signals  for  the  generation  of  the  RSNs 
(Brookes  et  al.,  2011;  Pasquale  et  al.  2010).  This  conforms  well  with  findings  that  low- 
frequency  oscillations  coordinate  long-range  communication  (Stein,  Chiang,  and  Konig, 
2000).  However,  these  MEG  findings  do  not  align  entirely  with  the  role  of  gamma 
oscillations  in  local  neural  activity  (Buzsaki  and  Wang,  2012)  and  as  counterparts  of  BOLD 
signaling  (Logothetis  et  al.,  2001). 

Cross-frequency  coupling,  in  particular  phase-amplitude  coupling,  has  been  attributed  an 
important  role  for  communication  between  brain  areas  (Canolty  et  al.,  2006;  Canolty  and 
Knight,  2010;  Jensen  and  Colgin,  2007).  Therefore  cross-frequency  coupling  between  a  low- 
frequency  phase  and  gamma  activity  may  be  an  ideal  candidate  for  the  communication  in  the 
RSN.  Based  on  these  previous  results  we  here  test  two  hypotheses.  First,  interregional 
correlated  fluctuations  during  rest  are  regulated  by  the  coupling  between  the  phase  of  slow- 
frequency  activity  and  the  amplitude  of  high-frequency  oscillations.  Second,  this  mechanism 
defines  the  principal  modes  of  resting-state  connectivity. 

To  test  these  hypotheses  we  use  a  model  (megPAC),  which  consists  of  the  interpolated 
gamma  amplitude  at  key-events  of  the  strongest  coupling  low-frequency  phase.  The 
rationale  for  this  signal  model  posits  that  the  excitability  cycles  of  local  populations  need  to 
be  synchronized  to  form  a  network.  Therefore,  a  core  hypothesis  is  that  the  phases  of  local 
low-frequency  oscillations  emerging  from  communicating  brain  regions  need  to  be  time- 
locked,  with  no  phase  delay.  This  condition  is  based  on  evidence  from  the  small-world 
architecture  of  brain  connectivity  suggesting  that  cortical  or  subcortical  hubs  facilitate 
synchronization  of  distant  oscillatory  activity  with  zero-lag  (i.e.  with  no  phase  delay) 

(Haider  et  al.,  2006).  Population  excitability  cycles  commonly  occur  at  lower  frequencies 
(Osipova,  Hermes,  and  Jensen,  2008).  Therefore  this  slower  oscillatory  rhythm  may  provide 
a  gating  mechanism  that  time-marks  the  operations  of  local  circuits,  revealed  by  bursts  of 
higher-frequency  activity  (Buzsaki  and  Draguhn,  2004;  Buzsaki  and  Wang,  2012;  Lakatos  et 
al.,  2005).  Along  this  hypothesis,  an  additional  condition  to  network  formation  requires 
coherent  high-frequency  oscillations  between  regions:  this  would  facilitate  local  post- 
synaptic  integration  and  spiking  activity  to  occur  in  concert  amongst  network  elements. 
Previous  observations  suggest  that  high-frequency  bursts  are  preferentially  occurring  about 
the  trough  of  the  low-frequency  phase  cycles  of  local  field  potentials  (Canolty  et  al.,  2006). 

Consequently,  we  confirm  and  extend  non-invasively  in  healthy  participants  for  the  entire 
brain  previous  observations  of  local  phase-amplitude  coupling  (PAC)  (Canolty  et  al.,  2006; 
Osipova,  Hermes,  and  Jensen,  2008).  We  show  that  the  oscillatory  components  conditioning 
the  preferred  timing  of  high-frequency  activity  through  PAC  span  a  relatively  wide  low- 
frequency  range:  from  delta  (2-4  Hz),  theta  (4-8  Hz),  to  alpha  (8-12  Hz)  bands.  With  these 
results  we  demonstrate  that  PAC  contributes  to  the  network  connectivity,  thereby  delivering 
a  principle  for  global  communication  across  the  brain  at  rest. 
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2  Materials  and  Methods 

Note  that  most  data  preprocessing  and  MEG  source  imaging  was  performed  using 
Brainstorm  (Tadel  et  al.,  2011).  Brainstorm  is  an  open-source  software,  freely  available  to 
the  academic  community  (http://neuroimage.usc.edu/brainstorm/).  All  implementation 
details  are  therefore  readily  documented  and  can  be  verified  in  Brainstorm's  code. 

Data  Acquisition 

The  study  was  approved  by  the  local  ethics  committee,  in  accordance  with  the  Declaration 
of  Helsinki.  12  healthy,  right-handed  subjects  (4  females,  8  males;  age  range:  21-41  y.o.) 
were  recruited  to  participate  in  the  study  and  all  subjects  gave  informed  consent. 

The  participants  were  tested  for  possible  magnetic  artifacts  in  a  short  preliminary  MEG  run. 
Subjects  preparation  consisted  of  taping  3  to  4  head-positioning  coils  on  the  subject's  scalp. 
The  positions  of  the  coils  were  measured  relative  to  the  subject's  head  using  a  3-D  digitizer 
system  (Polhemus  Isotrack).  To  facilitate  anatomical  registration  with  MRI,  about  100 
additional  scalp  points  were  also  digitized.  One  pair  of  electrodes  was  positioned  and  taped 
across  the  participants'  chest  (one  above  the  inferior  extremity  of  the  left  rib  cage  and  one 
over  the  right  clavicle)  to  capture  electrocardiographic  (ECG)  activity  synchronized  with  the 
MEG  traces.  Similarly,  one  pair  of  electrodes  was  attached  above  and  below  one  eye  to 
detect  eye-blinks  and  large  saccades  (EOG). 

5  subjects  were  measured  in  seated  position  with  the  Elekta-Neuromag  VectorView  306- 
channel  system  with  a  sampling  rate  of  2000  Hz  (0.03  Hz  high-pass  online  filter,  660  Hz 
anti-aliasing  low-pass  online  filter);  7  subject  were  measured  in  seated  position  using  the 
275-channel  VSM/CTF  system  with  a  sampling  rate  of  2400  Hz  (no  high-pass  filter,  660  Hz 
anti-aliasing  online  low-pass  filter).  Magnetic  shielding  was  provided  by  magnetically- 
shielded  rooms  (MSR)  with  full  3-layer  passive  shielding  for  the  CTF/VSM  system,  and 
single-layer  shielding  with  Maxfilter  active  flux-compensation  for  the  Elekta-Neuromag 
system.  The  combination  of  the  two  recording  systems  should  not  influence  our  results, 
because  all  calculations  and  combination  of  the  results  are  performed  at  the  source  level. 
Moreover  for  task  related  studies  it  was  shown  that  the  different  MEG  systems  yield 
essentially  identical  results  (Weisend  et  al.,  2007). 

At  the  beginning  of  each  MEG  run,  the  location  of  the  subject's  head  within  the  MEG 
helmet  was  measured  by  energizing  the  head-positioning  coils,  following  standard 
procedures.  For  each  subject,  between  5  and  30  minutes  (20  minutes  on  average  for  the 
subject  group)  of  MEG  data  were  acquired,  during  an  average  of  5  runs  of  2  to  10-minute 
duration.  The  only  instruction  given  to  the  participants  was  to  keep  their  eyes  open  and  to 
relax  without  falling  asleep. 

A  2-minute  empty-room  recording,  with  the  same  acquisition  parameters,  and  with  no 
subject  present  in  the  MSR,  was  used  to  capture  some  of  the  sensor  and  environmental  noise 
statistics,  which  were  used  in  the  source  estimation  process,  as  explained  below. 
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For  subsequent  cortically-constrained  MEG  source  analysis,  a  Tl-weighted  MRI  acquisition 
of  the  cerebrum  was  obtained  from  each  participant  either  at  least  one  month  before  the 
MEG  session  or  afterwards. 

Data  Pre-processing 

MEG  traces  were  pre-processed  to  verify  data  quality  and  to  reduce  contamination  from 
artifacts  (cardiac,  eye  movements  and  blinks,  environmental  noise).  Data  from  the  Elekta- 
Neuromag  system  were  preprocessed  using  signal-space-separation  (SSS)  (Taulu,  Kajola, 
and  Simola,  2004),  as  implemented  in  the  Maxfilter  noise  reduction  system  from  Elekta- 
Neuromag.  Default  SSS  settings  were  used:  orders  of  spherical  harmonic  expansions  for  the 
inner  and  outer  source  models  were  8  and  3,  respectively.  Data  from  the  CTF/VSM  system 
were  corrected  with  the  manufacturer's  3rd  order  gradient  compensation  system  (no 
parameter  setting  required).  The  projectors  obtained  were  propagated  to  the  corresponding 
MEG  source  imaging  operator. 

All  recordings  were  visually  inspected  to  detect  segments  contaminated  by  head  movements 
or  remaining  environmental  noise  sources,  which  were  discarded  from  subsequent  analysis. 
Heart  and  eye  movement/blink  contaminations  were  attenuated  by  designing  signal-space 
projections  (SSP)  from  selected  segments  of  data  about  each  artifactual  event  (Nolte  and 
Curio,  1999).  Using  Brainstorm's  ECG  and  EOG  detection  functionality  (Tadel  et  al.,  201 1). 
heartbeat  events  were  automatically  detected  at  the  R  peak  of  the  ECG's  QRS  complex,  and 
eye  blink  events  were  determined  automatically  at  the  peaks  of  the  EOG  traces.  Projectors 
were  defined  using  principal  component  analysis  (PCA)  of  these  data  segments  filtered 
between  10  and  40  Hz  (for  heartbeats)  or  1.5  and  15  Hz  (for  eye  blinks)  in  a  160-ms  time 
window  centered  about  the  heartbeat  event,  or  400  ms  around  the  eye  blink  event.  The 
principal  components  that  best  captured  the  artifact’s  sensor  topography  were  manually 
selected  as  the  dimension  against  which  the  data  was  orthogonally  projected  away  from,  also 
using  the  routines  available  in  Brainstorm.  Note  that  in  most  subjects,  the  first  principal 
component  was  sufficient  to  attenuate  artifact  contamination.  The  projectors  obtained  for 
each  subject  were  propagated  to  the  corresponding  MEG  source  imaging  operator  as 
explained  below.  Powerline  contamination  (main  and  harmonics)  was  reduced  by  complex 
match  filtering  with  1-Hz  resolution  bandwidth  for  sinusoidal  removal,  also  available  in 
Brainstorm.  The  pre-processed  data  were  resampled  at  1000  Hz,  using  the  polyphase  filter 
implementation  from  Matlab  (The  Mathworks,  MA,  USA)  with  default  parameters. 

The  scalp  and  cortical  surfaces  were  extracted  from  the  MRI  volume  data.  A  surface 
triangulation  was  obtained  for  each  envelope  using  the  segmentation  pipeline  available  in 
Brainvisa  (Riviere  et  al.,  2009)  (http://brainvisa.info),  with  default  parameter  settings  and 
subsequently  imported  into  Brainstorm.  The  individual  high-resolution  cortical  surfaces 
(about  75,000  vertices  per  surface)  were  down-sampled  to  about  15,000  triangle  vertices 
(also  with  a  Brainstorm  process)  to  serve  as  image  supports  for  MEG  source  imaging. 

MEG  Source  Imaging 

Forward  modeling  of  neural  magnetic  fields  was  performed  using  the  overlapping-sphere 
technique  implemented  in  Brainstorm  (Huang,  Mosher,  and  Leahy,  1999).  In  this  method. 
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one  sphere  is  automatically  adjusted  locally  to  the  individual  scalp  surface  under  each 
magnetic  sensor  to  compute  the  corresponding  lead  field  analytically.  This  method  has  been 
shown  to  provide  the  best  trade-off  between  modeling  precision  and  numerical  accuracy 
(Huang,  Mosher,  and  Leahy,  1999).  The  lead-fields  were  computed  from  elementary  current 
dipoles  distributed  perpendicularly  to  the  cortical  surface  from  each  individual  (Baillet, 
Mosher,  and  Leahy,  2001). 

MEG  source  imaging  was  performed  by  linearly  applying  the  weighted-minimum  norm 
operator  W  onto  the  preprocessed  data: 


W=At(AAt+AC)  1  (l) 


A  is  the  gain  matrix  from  the  forward  solution  and  C  the  spatial  covariance  of  the  noise  in 
the  recordings.  The  empirical  estimate  of  C  was  obtained  from  the  empty-room  recording 
described  above  (Baillet,  Mosher,  and  Leahy,  2001).  X  is  a  scalar  regularization  parameter. 
The  weights  of  the  imaging  operator  attenuate  the  decrease  of  signal-to-noise  ratio  between 
superficial  and  deeper  cortical  sources,  and  between  tangential  and  radially-oriented  sources 
with  respect  to  each  of  the  individual  spheres. 

As  mentioned  before,  the  data  were  projected  away  from  the  spatial  components  of  artifact 
contaminants.  For  consistency  between  the  projected  data  and  the  model  of  their  generation 
by  cortical  sources,  the  forward  operator  was  itself  projected  away  from  the  same 
contaminants  using  the  same  projector  as  for  the  MEG  data.  This  procedure  is  automatically 
verified  and  applied  in  Brainstorm. 

(-Amplitude  Coupling 

An  analytic  measure  of  cross-frequency  coupling  was  computed  from  each  of  the 
elementary  MEG  source  time  series,  in  each  individual  subject.  For  all  following 
computations,  all  artifact  free  recordings  from  each  subject  were  used.  They  were  therefore 
concatenated  in  the  source  space.  We  used  the  direct  phase-amplitude-coupling  (PAC) 
measure  (Ozkurt  and  Schnitzler,  2011),  defined  as  follow: 


(2) 


N  is  the  data  length,  n  a  running  time  sample;  a  is  the  signal  amplitude  at  the  high  frequency 
fa,  (p  is  the  signal  phase  at  the  low  frequency /^.  The  signal  phase  and  amplitude  were 
computed  using  the  chirplet  transform  (Mann  and  Haykin,  1991),  following  similar 
methodology  as  Canolty  et  al.  (2006). 

At  each  of  the  ~15,000  individual  source  locations  and  for  each  subject,  the  maximum  PAC 
score  was  computed  from  all  frequency  pairs  with G  [2,48  ]  Hz  and/a  G  [80,150] 
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Hz.  The  binning  width  for was  0.5  Hz  from  2  to  12  Hz  and  2  Hz  from  12  to  48  Hz.  The 
bins  of  fa  were  logarithmic allly  spaced.  The  optimal  low  frequency  phase  in  terms  of  largest 
PAC  score  with/,  €  [80,  150]  Hz  was  determined  at  each  cortical  location: 


(f<p,  /Q)=arg  max  PACify,  fa)-  (3) 

U,fa  () 


Statistical  analysis  of  Phase-Amplitude  Coupling 

Because  no  analytical  significance  measure  is  available  that  takes  into  account  the  usual 
temporal  correlation  in  neural  signals,  a  bootstrapping-based  approach  is  used  to  generate 
the  distribution  of  the  PAC  measure  under  the  null  hypothesis  HO  of  “no  PAC  and  serial 
correlation  of  the  type  found  in  neural  signals”.  In  particular,  we  generated  PAC  scores 
under  HO  using  the  following  procedure:  For  each  subject,  100  random  time  series  with  1  If 
characteristics  were  generated  with  durations  identical  to  the  subject's  actual  recordings'. 

The  l//- random  time  series  were  generated  by  applying  an  AR  filter  to  normally  distributed 
random  time  series  (Kasdin,  1995).  To  evaluate  whether  the  pink  noise  data  is  really  able  to 
capture  the  salient  autocorrelation  structure  of  the  data,  we  compared  autocorrelograms  of 
the  amplitude  in  different  frequency  ranges  in  the  experimental  and  noise  data  and  found 
them  to  be  quite  close  to  each  other  (summarized  e.g.  by  the  half-life  of  the  correlation). 

Based  on  these  artificial  noise  time  series,  the  critical  values  for  significance  testing  were 
computed  the  following  way:  for  each  random  time  series,  the  PAC  value  obtained  for  each 
low  frequency  phase and  high  frequency  amplitude  fa  £  [80,  150]  Hz  was  calculated  in 
the  same  manner  as  for  the  experimental  data  times  series.  From  the  100  surrogate 
repetitions,  the  value  at  the  95th  quantile  of  PAC  scores  for  every  low/high  frequency  pair 
was  determined  as  the  significance  threshold  at  p  <  0.05  for  actual  recordings  with 
Bonferroni  correction  (15,000  sources  and  429  frequency  pairs,  using  linear  interpolation 
when  necessary). 

The  optimal  low  frequency  for  phase/' ^  was  classified  in  one  of  the  typical  sub-bands  of 
electrophysiological  neural  signals:  S=  [2,4]  Hz,  6=  [4,8]  Hz,  a  =  [8,12]  Hz,  j3=  [12,30]  Hz 
or  y=  [30,48]  Hz.  Every  f' q  value,  obtained  at  all  cortical  locations,  was  mapped  onto  the 
individual  cortical  surface  from  each  subject.  For  group  analysis,  the  resulting  individual 
significant  PAC  maps  were  projected  onto  the  Colin27  brain  template,  using  Brainstorm's 
multilinear  registration  procedure.  The  resulting  anatomically-registered  maps  were  then 
averaged  across  all  12  subjects  excluding  the  non-significant  regions  and  displayed  in  Fig. 
2a. 

Estimation  of  frequency  with  maximal  power 

For  every  subject  at  every  cortical  location  we  calculated  the  power  spectral  density  using 
the  fast  fourier  transform  with  half  overlapping  4096  data-points  long  Hanning  windows. 
This  yielded  a  frequency  resolution  of  0.24Hz.  We  took  the  maximal  frequency  from 
2-48Hz  in  accordance  to  the  maximal  PAC  determination  and  calculated  the  difference  to 
the  maximal  PAC  frequency  for  all  sources. 
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Generation  of  Figure  2b 

The  method  used  to  obtain  Figure  2b  was  akin  to  that  of  Canolty  et  al.,  2006  and  applied  to 
the  source  signals  at  a  random  selection  of  cortical  locations.  First,  the  troughs  at  the  optimal 
nesting  frequency  /' ^  were  determined  from  the  signal  phase  time  series  (pn(f  E  {-tt,  7r]  - 
with  (pn(f*(p )  =  7r  corresponding  to  a  trough  -  obtained  with  chirplets  (Mann  and  Haykin, 
1991).  Second,  the  amplitude  a(n,f)  of  the  source  signals  with/E  [40,  250]  Hz  was 
determined,  also  using  the  chirplet  transform.  Third,  the  amplitude  time  series  a(n,  f)  was  z- 
scored  for  each  frequency  bin: 


a(n,f) 


a(n,  f)  -  { a{n,f )) 

“(/) 


(4) 


where  ( a(f ))  is  the  time  average  of  a{n,j),  and  d(f)  is  the  temporal  standard  deviation  of 
a(n,f),  estimated  over  all  N  time  samples,  for  each  subject.  Fourth,  the  time-frequency 
decomposition  d(n,f)  was  epoched  over  the  time  window  of  [-0.5,  0.5]  s  about  each  of  the 
previously  identified  phase  troughs  yielding  d(k,  n,f).  The  number  of  epochs  k  ranged 
between  2,600  and  24,000  depending  on  the  data  length  for  each  subject.  Note  that  due  to 
the  data  concatenation  of  the  different  runs  for  each  subject,  there  are  slight  jumps  in  the 
data.  However,  these  occur  at  most  10  times  per  subject  so  that  their  influence  should  be 
negligible.  Fifth,  time-frequency  representations  of  a(k,  n,  f)  were  averaged  across  all  k 
epochs  as  shown  in  Figure  2b.  Sixth,  the  source  signals  were  also  averaged  across  all  k 
epochs  to  obtain  (a(k,  n ))  the  time  series  shown  at  the  bottom  of  the  time-frequency  maps  in 
Figure  2b. 

The  significance  of  phase-amplitude  coupling  in  terms  of  correlated  temporal  modulations 
between  the  high-frequency  nested  components  and  the  low-frequency  nesting  oscillations 
over  the  epoch  time  window  was  determined  by  evaluating  the  slope  between  the  event- 
related  average  of  the  broadband  cortical  time  series,  conditioned  to  phase  trough  events 
(<Pn(Ttf>=  /)>  and  the  average  of  the  time  variations  d(k,  n,  f)  at  every  frequency  bin/.  The 
slope  coefficient  corresponds  to  the  correlation  coefficient.  To  estimate  the  slope  we  used  an 
ordinary  least  square  estimate.  Using  Newey  and  West,  1987  robust  standard  errors  and  the 
asymptotically  normal  distribution  of  the  coefficients,  we  assessed  the  significance  of  the 
slope  coefficients.  Significance  thresholds  were  set  at  each  frequency  bin/at  p  <  0.05,  with 
Bonferroni  correction. 


megPAC  and  control  models 

From  the  MEG  source  signals  we  generated  a  signal  model  to  evaluate  factors  of  long-range 
connectivity  between  brain  regions  in  the  resting-state.  The  generation  of  the  megPAC 
signal  models  is  illustrated  schematically  in  Figure  1 :  the  gamma  amplitude  of  the  source 
signal  in  the  [80,  150]  Hz  was  linearly  interpolated  with  Malab's  interpl  function  to  the  time 
occurrences  of  the  peaks  and  troughs  of  the  optimal  low-frequency  phase  //>  yielding  the 
largest  PAC  within  the  [80,150]  Hz  high-gamma  range,  as  obtained  from  Eq.  3.  The  gamma 
amplitude  in  the  [80,  150]  Hz  range  was  obtained  from  the  chirplet  transform  of  the  original 
source  signal,  as  explained  above.  To  compensate  for  the  \/f  decrease  in  signal  power,  the 
gamma  amplitude  was  calculated  at  logarithmically-spaced  frequency  bins  and  each 


Neuroimage.  Author  manuscript;  available  in  PMC  2016  May  01. 


Author  Manuscript  Author  Manuscript  Author  Manuscript  Author  Manuscript 


Florin  and  Baillet 


Page  8 


temporal  variation  at  each  bin  was  z-transformed.  Therefore  the  values  of  the  mean  and 
standard  deviation  were  taken  over  the  entire  individual  recording  (N  time  samples)  for  each 
frequency  bin.  An  estimation  of  the  signal  amplitude  in  the  [80,  150]  Hz  range  was  obtained 
from  averaging  over  all  frequency  bins  in  that  range.  The  resulting  megPAC  (tf)  values 
were  subsequently  downsampled  to  10  Hz  after  applying  an  anti-aliasing  filter. 

To  verify  the  validity  of  the  megPAC  model,  we  generated  three  alternative  models: 

1.  We  investigated  whether  identifying  the  optimal  coupling  frequency/' ^  at  each 
cortical  location  was  critical  to  the  extraction  of  resting-state  networks.  We 
therefore  generated  another  signal  model  following  the  same  approach  as  with  the 
megPAC  model,  but  with  a  nesting  frequency  arbitrarily  set  in  the  Grange  at  5.4 
Hz,  for  each  of  the  15,000  tested  brain  locations,  as  used  by  Canolty  et  ah,  2006. 

2.  We  tested  whether  the  identification  of  certain  low-frequency  events  is  essential  at 
all.  Therefore  we  selected  random  time-points  of  the  gamma  amplitude  and 
interpolated  between  those. 

3.  We  generated  12  sets  of  random  time  series  with  Ilf  characteristics,  which 
simulated  MEG  sensor  activity.  These  data  sets  were  source  reconstructed  with  a 
weighted  minimum-norm  operator  on  the  Colin27  brain.  For  these  data  the  optimal 
PAC  frequency  and  the  megPAC  model  were  calculated.  This  model  was  used  to 
determine  effects  in  the  obtained  networks,  from  the  processing  pipeline  and  the 
sensitivity  profile  of  MEG  source  imaging. 

Extraction  of  Resting-State  Networks  and  Statistical  Thresholds 

The  following  procedure  for  extracting  resting-state  networks  as  dominant  patterns  of 
connectivity  between  cortical  sources  was  applied  to  each  signal  model  introduced  in  the 
previous  section: 

1.  Projection  onto  anatomical  template:  The  individual  megPAC  time  series  were 
projected  onto  the  Colin27  cortical  surface  template  (consisting  of  N  ,j=  1 5,000 
triangle  vertices)  using  the  multilinear  scaling  and  registration  procedure. 

2.  Spatial  smoothing  and  signal  concatenation:  The  cortical  maps  obtained  from  the 
megPAC  and  alternative  signal  models  were  spatially  smoothed  with  a  non¬ 
recursive  kernel  of  7mm  full-width  at  half  maximum  applied  at  the  surface  of  the 
cortex.  The  time  series  of  all  subjects  were  concatenated,  following  previous 
reports  in  MEG  (Brookes  et  al.,  201 1)  and  fMRI  (Smith  et  ah,  2009)  RSN  analysis. 

3.  Connectivity  matrices:  The  pairwise  correlation  coefficients  were  calculated 
between  the  time  series  taken  at  each  pair  of  cortical  sources,  yielding  the  Nv  x  Nv 
array  C  of  empirical  correlation  coefficients. 

4.  Dimensionality  reduction:  We  conducted  the  extraction  of  RSN  maps  by  obtaining 
the  principal  spatial  modes  of  C.  To  make  this  operation  computationally  tractable, 
each  row  of  C  was  projected  orthogonally  onto  those  of  C1175,  a  subset  of  1 175 
rows  from  C  corresponding  to  cortical  locations  evenly  distributed  over  the 
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template  cortical  surface.  This  approach  was  proposed  in  RSN  MRI  by  Yeo  et  al., 
2011  and  yielded  the  Nv  x  1 175  reduced  connectivity  array  P,  defined  as  follows: 

P— CC1175,  (5) 


where  T  denotes  matrix  transpose. 

5.  Compensation  for  MEG-imaging  resolution  and  sensitivity:  To  minimize  the 
influence  of  the  uneven  sensitivity  of  MEG  source  imaging  to  anatomical  features 
(source  orientation  and  depth),  the  previous  steps  were  applied  to  surrogate  MEG 
sensor  data  (independent  and  identically  distributed,  zero-mean,  unit  variance 
Gaussian  sensor  time  series  of  same  length  as  the  original  data's,  for  12  surrogate 
individual  data).  The  corresponding  pairwise  reduced  connectivity  array  for 
surrogate  data  was  obtained  and  denoted  P.  From  the  singular  value  decomposition 
(S  VD)  of  P,  P  =  U  S  VT,  the  first  spatial  singular  mode  in  U  was  used  to  create  II, 
the  orthogonal  projector  away  from  the  principal  mode  of  spurious  connectivity  due 
to  MEG-imaging  resolution  and  uneven  sensitivity. 

6.  RSN  maps  as  principal  modes  of  connectivity:  The  RSN  maps  were  defined  as  the 
principal  spatial  modes  of  P  after  they  were  corrected  using  the  PII  projector  (see 
previous  step),  and  identified  as  U,  obtained  from  the  following  SVD: 

Pf[=USVT.  (6) 

Figure  3  shows  the  cortical  maps  thresholded  at  40%  of  their  maximum  value  and  a 
minimum  2.5  cm2  cluster  size,  of  5  principal  spatial  modes  of  the  signal  CSM  projected 
away  from  the  principal  noise  CSM  pattern. 

The  extracted  network  across  all  subjects  of  the  megPAC  and  the  alternative  models,  which 
corresponded  best  to  the  DMN,  were  compared  in  Figure  4b.  We  used  the  default-mode  map 
as  shown  in  the  lower  right  part  of  Figure  3  as  documented  in  Vincent  et  al.,  2008  to  test  for 
the  correspondence  between  the  DMN  map  and  the  maps  based  on  the  models.  We 
calculated  the  ratio  of  correctly  detected  DMN  areas  versus  false  positively  detected  areas. 
The  ratios  for  the  models  were  compared  with  students  t-test  and  subsequent  Bonferroni 
correction.  The  ratio  allows  for  estimating  the  correspondence  in  the  spatial  domain. 

For  each  seed  location  /,  the  z-scored  connectivity  pattern  of  the  seed  region  (i.e.,  the  /th  row 
of  C,  standardized  with  respect  to  the  mean  and  standard  deviations  observed  with  surrogate 
time  series  at  the  same  location  /)  was  obtained.  All  other  brain  locations  were  labeled 
according  to  the  maximum  observed  connectivity  z-score  amongst  the  seed  candidate 
region.  Only  regions  with  z-scores  above  7  were  considered  for  labeling.  This  method  was 
used  to  generate  Figure  4a,  where  seeds  of  interest  were  positioned  in  3  distinct  longitudinal 
aspects  of  the  inferior  frontal  gyrus. 

All  maps  were  spatially-interpolated  onto  a  higher-resolution  version  of  the  template  brain, 
with  120,000  vertices,  using  the  procedure  available  in  Brainstorm. 


Neuroimage.  Author  manuscript;  available  in  PMC  2016  May  01. 


Author  Manuscript  Author  Manuscript  Author  Manuscript  Author  Manuscript 


Florin  and  Baillet 


Page  10 


3  Results 

Ubiquitous  Phase-Amplitude  Coupling  across  the  Human  Cortex 

We  first  tested  whether  physiologically  relevant  PAC  can  be  detected  globally  from  human 
brain  dynamics  during  rest.  To  this  end,  we  performed  source  imaging  (Baillet,  Mosher,  and 
Leahy,  2001;  Tadel  et  al.,  2011)  from  the  MEG  recordings.  At  every  cortical  location,  the 
direct  PAC  measure  (Ozkurt  and  Schnitzler,  2011)  was  computed  between  the  low- 
frequency  phase  and  high-frequency  amplitude  of  oscillatory  MEG  source  signals.  Multiple 
candidate  low-frequency  components  for  phase  were  tested  between  2  to  12  Hz  in  bins  of 
0.5  Hz  and  from  12  to  48  Hz  in  2  Hz  steps.  The  reason  for  testing  systematically  a  wide 
range  of  low-frequency  components  is  the  lack  of  consensus  across  previous  published 
findings  about  a  preferred  band  of  low-frequency  oscillations  for  phase  coupling  with  high- 
gamma  amplitude  bursts.  Phase-locked  high-frequency  amplitude  fluctuations  were  tested  in 
the  higher  gamma  frequency  range  (80-150  Hz).  For  every  source,  direct  PAC  scores  were 
obtained  for  each  tested  pair  of  low-frequency  phase  and  high-frequency  amplitude.  We  then 
determined  the  frequency  pair  with  the  highest  direct  PAC  score  and  assessed  statistical 
significance. 

For  the  highest  PAC  score  reported  here  41%-61%  of  the  sources  at  the  individual  level 
showed  significant  PAC  (p  <  0.05,  Methods).  The  main  modes  across  subjects  of  the  low- 
frequency  phase  components  were  prominently  in  the  delta  (2-4  Hz)  (65.3%  of  the  cortical 
surface  area  averaged  across  subjects,  with  18.6-34.8%  HR)  and  theta  (4-8  Hz)  bands 
(33.7%  of  the  cortical  surface  area  averaged  across  subjects,  with  14.3-20.0%  HR).  The 
alpha  (8-12  Hz)  band  was  found  as  low-frequency  phase  component  over  1%  (3. 9-7. 8% 

HR)  of  the  cortex.  Low-frequency  phase  components  in  the  beta  band  (12-30  Hz)  remained 
spatially  marginal  (0%  on  average,  with  0.1-10.3%  HR).  No  significant  PAC  with  low- 
frequency  phase  in  the  low-gamma  range  (30-48  Hz)  was  found.  Figure  2a  shows  the  low 
frequency  phase  of  the  PAC  averaged  across  subjects.  For  this  average  the  low  frequency 
phase  of  the  maximal  and  significant  PAC  was  used.  Note  that  no  spatial  pattern  was  found 
in  these  results.  It  rather  demonstrates  that  delta  and  theta  are  the  main  low-frequency 
phases.  Having  identified  PAC  over  the  whole  brain  suggests  that  it  is  not  confined  to 
certain  brain  regions,  but  may  represent  a  fundamental  building  block  for  communication 
within  the  brain. 

To  ensure  that  the  maximal  PAC  is  not  driven  by  the  power  of  the  low-frequency  phase  we 
determined  the  difference  in  frequency  of  the  low-frequency  phase  from  the  maximal  PAC 
and  the  frequency  corresponding  to  the  maximal  power.  In  80%  of  all  sources  across 
subjects  these  two  frequencies  were  at  least  1Hz  apart  and  in  69%  at  least  2Hz. 

To  provide  an  illustrative  example  of  the  phase-amplitude  coupling,  we  used  an  event- 
related  average  of  spectrograms  of  ongoing  cortical  signals  obtained  about  the  troughs  of  the 
dominant  low-frequency  phase  component.  This  approach  is  similar  to  the  one  of  Canolty  et 
al.,  2006.  These  results  are  shown  for  arbitrarily  chosen  brain  locations  of  one  subject  in 
Figure  2b.  We  found  significant  levels  of  temporal  correlation  (p  <  0.05,  Methods)  between 
the  amplitude  fluctuations  of  high-gamma  activity  and  the  phase  cycles  of  local  low- 
frequency  oscillations. 
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Phase-Amplitude  Coupling  and  Resting-State  Networks 

We  tested  the  second  hypothesis  that  the  ongoing  dynamics  of  local  PAC  oscillations 
determine  network  formation  between  the  nodes  of  RSN. 

In  order  to  capture  the  most  salient  features  of  ongoing  local  PAC  dynamics  from  ongoing 
MEG  source  time  series,  a  signal  transform  was  applied  at  each  cortical  location:  The 
amplitude  of  high-gamma  (80-150  Hz)  bursts  was  linearly  interpolated  between  the  troughs 
and  peaks  of  the  low-frequency  phase  cycles  with  highest  PAC  to  high-frequency  amplitude. 
In  essence,  the  interpolation  of  the  peaks  and  troughs  reduces  the  original  MEG  source  time 
series  to  the  peaks  and  troughs  of  optimal  PAC  cycles,  which  are  key  timing  events  for  long- 
range  communication.  In  local  field  potential  recordings,  the  convention  is  that  these  key 
events  occur  at  the  troughs  of  the  low-frequency  phase  (Buzsaki  and  Draguhn,  2004; 

Canolty  et  ah,  2006).  However,  MEG  source  estimation  detects  current  flows,  not  electrical 
potentials,  and  presents  a  +  /r  phase  ambiguity,  which  depends  on  conventional  default 
orientation  of  elementary  current  dipole  sources  constrained  to  the  local  cortical  anatomy 
(Baillet,  Mosher,  and  Leahy,  2001).  Hence  we  considered  both  peaks  and  troughs  as  phase 
events  of  interest.  We  called  the  resulting  sub-sampled  signal  model  “megPAC”. 

megPAC  signals  were  obtained  for  each  individual  subject  and  processed  for  RSN 
extraction  (Methods).  Figure  3  shows  that  the  principal  spatial  modes  of  resting-state 
connectivity  obtained  with  the  megPAC  signal  model  show  similarities  to  the  typical  RSN 
reported  by  fMRI  studies. 

Our  results  therefore  confirm  that  RSN  are  determined  by  correlated  dynamics  of  local  PAC 
fluctuations  between  regional  nodes. 

We  further  tested  whether  more  fine-grained  and  anatomy-specific  analysis  of  functional 
connections  in  the  brain  can  be  detected  from  the  megPAC  model  with  a  seed-based 
analysis.  Anatomical  seeds  were  placed  in  the  orbitalis  (IFGor),  triangularis  (IFGtr),  and 
opercularis  (IFGop)  aspects  of  the  left  inferior  frontal  gyrus.  The  functional  connectivity 
patterns  obtained  were  strikingly  similar  to  those  reported  from  diffusion-weighted  MRI 
tractography  studies  and  post-mortem  dissections  (Catani  et  al.,  2012)  (Figure  4a): 
Connections  with  IFGtr  were  circumscribed  about  the  middle  lateral  frontal  and  central 
regions.  Connections  with  IFGop  extended  to  the  pre-supplementary  motor  area  (SMA), 
anterior  and  middle  cingulate  regions,  the  insula,  the  mouth  sensori-motor  cortex,  and  the 
superior  and  middle  temporal  lobe. 

Validation  of  the  signal  models 

We  quantified  the  correspondence  of  the  resting-state  networks  identified  from  using  the 
megPAC  model  and  the  three  alternative  models  described  in  the  methods  with  a  natural 
benchmark,  using  a  standard  DMN  map  from  fMRI  (Vincent  et  al.,  2008).  We  calculated  the 
ratio  of  the  percentage  correspondence  versus  percentage  falsely  identified  DMN  regions  for 
each  signal  model.  Only  the  megPAC  model  shows  a  ratio  larger  than  1,  indicating  that  it 
correctly  identified  a  higher  percentage  of  actual  DMN  regions  as  belonging  to  the  DMN 
than  non-DMN  regions.  Therefore  the  megPAC  model  yielded  a  significantly  higher 
correspondence  to  the  DMN  than  all  other  tested  models  (Figure  4b,  p  <  0.001).  Furthermore 
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the  megPAC  model  with  theta  only  and  the  megPAC  model  with  random  troughs  yielded  a 
lower  correspondence  for  the  DMN  than  the  megPAC  model  from  random  data.  This  clearly 
indicates  that  only  the  megPAC  model  properly  captures  the  dynamics  of  the  resting  state. 
For  the  DMN,  medial  prefrontal  regions,  the  anterior  temporal  lobe,  posterior  cingulate 
cortex  (PCC),  and  the  left  and  right  inferior  parietal  lobules  are  the  main  components  found 
with  resting-state  fMRI.  Our  analysis  essentially  revealed  the  prefrontal  aspects  of  the 
DMN,  with  no  identified  contribution  from  the  PCC  and  inferior  parietal  regions.  Note  that 
the  PCC's  contribution  to  the  DMN  was  also  not  detected  in  previous  MEG  studies  of  the 
RSN  using  data-driven  techniques  (Brookes  et  al.,  2011). 

We  calculated  the  ratio  of  corresponding  to  non-corresponding  identified  regions  for  the 
DAN,  the  visual  network,  and  the  right  fronto-parietal  network.  We  used  the  standard  maps 
shown  in  figure  3  (Vincent  et  ah,  2008;  Yeo  et  ah,  2011)  and  found  that  this  ratio  was  1.88 
for  the  DAN,  0.03  for  the  visual  network,  and  3.38  for  the  right  fronto-parietal  network.  The 
low  correspondence  for  the  visual  network  might  be  due  to  the  smaller  spatial  extension  of 
the  MEG-based  network.  The  MEG  analysis  detected  an  additional  medial  frontal 
component,  which  is  not  a  typical  part  of  the  fMRI-based  visual  network.  A  combination  of 
the  DAN  and  the  sensorimotor  network  (SMN)  was  detected  with  MEG.  This  is  in 
accordance  with  seed-based  analysis  in  fMRI  data  (Vincent  et  ah,  2008)  where  the 
DAN/SMN  is  also  extracted  when  placing  a  seed  in  the  superior  parietal  lobule  and  the 
middle  temporal  area  (MT+).  For  the  right  fronto-parietal  network,  the  same  lateral 
components  as  with  fMRI  were  detected  (prefrontal  cortex,  inferior  parietal  lobule,  and 
inferior  temporal  gyrus).  Yeo  et  ah  (2011)  described  additional  medial  components  for  the 
fronto-parietal  network,  which  were  not  detected  with  MEG.  Concerning  the  auditory 
network,  a  clear  activation  of  the  left  auditory  cortex  as  well  as  the  inferior  frontal  gyrus  was 
detected  with  MEG,  while  the  activation  of  the  right  auditory  cortex  was  missing. 

4  Discussion 

In  the  present  study  we  demonstrated  that  phase-amplitude  coupling  provides  a  mechanism 
for  brain  network  formation,  which  reconciles  previous  findings  and  theories  on  long-range 
communication  between  neural  populations. 

Our  first  result  is  that  local  PAC  can  be  demonstrated  during  rest  over  the  whole  cortex  and 
that  it  is  accessible  non-invasively  in  humans.  In  particular,  we  found  cortical  PAC  of 
ongoing,  local  activity  between  a  slower  low-frequency  phase  (predominantly  in  the  2-12Hz 
range)  and  bursts  of  high-frequency  activity  (80Hz  and  above)  at  the  peaks  and/or  the 
troughs  of  the  low-frequency  oscillations.  These  results  obtained  with  MEG  source  imaging 
extend  previous  reports  obtained  with  invasive  recording  techniques  and  restricted  to  a 
limited  number  of  brain  regions  (Canolty  et  al.,  2006;  Tort  et  al.,  2008,  2009;  Voytek  et  al., 
2010).  Thus,  a  comprehensive  map  of  PAC  measures  over  the  entire  cortical  surface  was 
obtained  in  each  tested  healthy  subject.  These  results  provide  a  promising  avenue  for 
systematically  studying  PAC  and  its  relevance  for  neural  communication  in  other  conditions 
than  rest. 
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Second,  we  showed  that  PAC  mechanisms  contribute  to  long-range  neural  communication 
between  brain  regions  in  the  resting-state.  For  this  purpose,  we  identified  the  principal 
modes  of  connectivity  from  the  correlated  fluctuations  of  a  signal  model  based  on  local 
dynamics  of  phase-amplitude  coupling  (megPAC).  Through  the  specific  sampling  of  the 
gamma-amplitude  with  the  restriction  to  the  peaks  and  troughs  of  the  low  frequency,  it  was 
ensured  that  both  features  of  the  PAC  were  represented  in  the  model.  Furthermore  only  if 
two  brain  regions  have  the  same  low -frequency-phase  and  gamma  amplitude,  will  the 
correlation  be  high.  Thereby  the  connectivity  maps  reveal  brain  regions  with  the  same  PAC 
characteristics.  The  megPAC  model  revealed  the  expected  major  resting-state  networks  and 
was  corroborated  with  a  seed-based  approach.  The  correspondence  to  a  standard  default  map 
was  demonstrated  through  statistical  comparison  of  corresponding  and  non-corresponding 
detected  regions.  Most  of  the  missing  regions  are  within  the  medial  aspects  of  the  cortical 
hemispheres,  which  might  be  due  to  the  lower  sensitivity  of  MEG  within  these  regions.  We 
do  not  have  the  individual  resting  state  fMRI  maps  and  therefore  a  direct  comparison  was 
not  possible. 

This  result  shows  that  it  is  possible  to  identify  motifs  of  megPAC  signal  correlation  that  are 
compatible  with  known,  finer  anatomical  connections  (Catani  et  al.,  2012).  Taken  together, 
our  results  suggest  that  phase-amplitude  coupling  is  an  important  component  of  the 
electrophysiology  underlying  network  generation  in  the  resting-state 

Furthermore,  our  results  reconcile  observations  from  previous  resting-state  MEG  studies  in 
humans  (Brookes  et  al.,  2011;  Pasquale  et  al.,  2010)  and  combined  BOLD  and  electro¬ 
physiology  recordings  in  animals  (Logothetis  et  al.,  2001).  They  demonstrate  that  both  the 
lower-frequency  components  reported  by  these  MEG  studies  and  the  high-frequency 
oscillations  found  in  invasive  recordings  are  actually  coupled  and  that  they  both  contribute 
to  a  unifying  mechanism  of  resting-state  brain  connectivity. 

The  fact  that  resting-state  networks  are  accessible  non-invasively  using  MEG  opens  new  and 
substantial  methodological  possibilities  for  investigating  the  dynamics  of  resting-state 
activity  at  multiple  temporal  and  spatial  scales.  In  particular,  using  the  time  resolution  of 
MEG  and  our  approach  for  the  analysis  of  resting-state  activity,  it  is  in  principle  now 
possible  to  study  the  dynamical  interplay  between  resting-state  networks  occurring  at  a  finer 
time  scale.  We  anticipate  these  results  will  stimulate  further  research  on  the  physiology  of 
the  cerebral  resting-state,  its  transition  in  response  to  stimuli  and  during  task  performance, 
and  will  ultimately  contribute  to  finding  new  markers  of  healthy  and  diseased  brain  function. 

Neural  Network  Formation  with  Synchronized  Gating 

Based  on  our  findings  here  and  previous  studies,  we  propose  the  model  of  synchronized 
gating  (SG).  SG  extends  two  previously  proposed  hypotheses  for  large-scale  neural 
communication:  the  communication-through-coherence  (CTC)  hypothesis  (Fries,  2005)  and 
the  binding-by-synchronization  (BBS)  model  (Singer,  1999;  Varela  et  al.,  2001).  The  idea  of 
the  synchronized  gating  model  is  that  the  synchronized  phase  of  slower  neural  oscillations 
provides  a  common  gating  mechanism  between  multiple  brain  regions:  the  low-frequency 
phase  time-marks  the  occurrence  of  coherent,  local  high-frequency  activity,  which 
comprises  regional  post-synaptic  integration  and  spiking  activity.  Underlying  this  model  is 
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the  previously  documented  finding  that  slower  oscillatory  rhythms  time-mark  the  operations 
of  local  circuits  (Buzsaki  and  Wang,  2012).  In  addition,  low-frequency  oscillations  have 
been  reported  to  coordinate  long-range  communication,  while  gamma  oscillations  (40-150 
Hz)  coordinate  local  communication  (Buzsaki  and  Draguhn,  2004;  Stein,  Chiang,  and 
Konig,  2000).  If  low-frequency  phases  are  asynchronous  (e.g.,  occurring  at  different 
frequencies  and/or  with  a  phase  delay),  the  “gate  is  closed”,  coherent  high-frequency 
oscillations  and/or  synchronous  cell  firing  are  less  likely  to  occur,  and  communication/ 
coactivation  between  two  brain  regions  is  suboptimal. 

Experimental  evidence  supported  by  computational  models  (Vicente  et  ah,  2008) 
demonstrates  that  precise  synchronization  of  neural  oscillations  can  emerge  with  zero  time 
lag  to  form  networks.  This  is  achieved  through  cortical  (Buzsaki  and  Wang,  2012;  Hipp  et 
ah,  2012)  or  subcortical  (Sherman  and  Guillery,  2002)  hubs,  which  act  as  dynamical  relays 
of  long-distance  cortico-cortical  communication.  This  network  architecture  circumvents  the 
influence  of  variable  axonal  lengths  and  conduction  delays  on  time  lags  between  neural 
signals  (Vicente  et  ah,  2008).  Thereby,  according  to  SG,  the  recruitment  of  regions  to 
participate  in  networks  becomes  a  dynamic  and  flexible  process,  because  the  recruitment  is 
not  conditioned  to  direct  connections  between  network  nodes.  The  flexibility  in  network 
formation  proposed  by  the  SG  model  could  therefore  be  essential  for  behavioral  adaptation 
and  promote  developmental  and  learning  abilities. 

In  summary,  network  formation  by  SG  requires  that  nodal  regions  demonstrate  synchronized 
low-frequency  activity  with  no  phase  delay,  which  regulates  the  occurrence  of  coherent 
levels  of  local  high-frequency  bursts.  In  essence,  it  encapsulates  the  essence  of  both  the  CTC 
and  BBS  hypotheses  in  a  single  proposition.  Figure  5  illustrates  the  basic  mechanisms 
posited  by  the  SG  model. 

Methodological  Considerations 

We  have  used  the  phase-amplitude  coupling  measure,  which  has  been  proposed  before  and 
yielded  the  most  robust  results.  Still,  there  are  known  methodological  pitfalls  with  any  PAC 
measure.  The  first  one  relates  to  the  phase-ambiguity  in  MEG  signals.  In  particular,  it 
remains  to  be  clarified  whether  the  high  frequency  bursts  occur  preferentially  at  the  trough, 
as  suggested  by  Canolty  et  ah,  2006,  or  another  phase  of  the  low-frequency  cycle. 

The  second  methodological  pitfall  concerns  the  sensitivity  to  low-frequency  phase 
components  with  asymmetric  cycles,  as  reported  by  local  recordings  (Buzsaki  and  Draguhn, 
2004:  Haider  et  ah,  2006;  Steriade  et  ah,  1996).  Related  to  this  edges  in  the  data  lead  to 
artificial  PAC  (Kramer,  Tort,  and  Kopell,  2008). 

Still  PAC  has  been  found  in  numerous  studies  using  ECoG,  LFP,  and  MEG  data  (Canolty  et 
ah,  2006;  Ozkurt  et  ah,  2011;  Tort  et  ah,  2008)  making  it  unlikely  that  PAC  is  artefactual. 
Moreover,  our  particular  exemplary  maps  in  Figure  2b  show  the  averaged  low-frequency 
signal.  In  this  average  no  asymmetry  of  the  oscillation  is  seen,  therefore  indicating  that  our 
results  are  not  driven  by  spurious  phase  asymmetry.  Furthermore  our  choice  for  chirplet 
time -frequency  decompositions  yielded  a  good  combination  for  time  and  frequency 
resolution. 
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The  third  methodological  issue  related  to  PAC  is  that  it  could  potentially  be  driven  by  the 
low-frequency  with  the  maximal  power.  Our  results  indicate  that  this  is  not  the  main  driving 
cause  of  the  PAC,  because  in  80%  of  all  subjects  these  two  frequencies  did  not  align  with 
one  another.  This  also  suggests  a  possible  alternative  mechanism  for  the 
electrophysiological  origins  of  resting-state  networks  than  proposed  in  previous  MEG 
studies  (Brookes  et  ah,  2011;  Hipp  et  ah,  2012;  Pasquale  et  ah,  2010). 

In  summary,  we  established  that  local  phase-amplitude  coupling  between  neural  oscillations 
is  a  fundamental  mechanism  for  long-range  communication  between  brain  regions.  With  a 
signal  model  based  on  phase-amplitude  coupling,  resting-state  networks  could  be  extracted. 
Based  on  these  results,  we  suggested  the  synchronized  gating  hypothesis  as  a  mechanism 
enabling  long-range  neural  communication  and  functional  connectivity  in  the  brain,  during 
rest  and  task  performance. 


A  Chirplet  transform 


Chirplets  separate  a  time  series  not  only  in  its  time  and  frequency  components,  but  also 
allow  to  consider  a  shear  in  time  and  frequency  as  well  as  a  time  dilation/frequency 
contraction.  To  construct  a  chirplet  a  windowing  function  such  as  a  Gaussian  window  is 
applied  to  a  chirp.  Therefore  in  the  case  of  a  Gaussian  windowing  function,  which  has  been 
used  in  the  present  paper,  the  chirplet  is  described  as  follows: 


(7) 


tc  is  the  time  center,  fc  the  frequency  center,  Ar  the  duration,  and  c  the  chirprate  of  the 
chirplet.  In  the  present  paper  the  chirprate  was  set  to  0.  These  parameters  give  a  greater 
flexibility  in  covering  the  time-frequency  domain.  For  an  in  depth  review  see  Mann  and 
Haykin,  1995. 
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Highlights 

•  Propose  and  experimentally  verify  synchronized  gating  as  mechanism  for 
communication. 

•  Non-invasive  detection  across  the  entire  cortex  of  local  phase-amplitude 
coupling. 

•  Spatial  modes  of  cross-frequency  coupling  correspond  to  resting-state  networks. 
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Figure  1.  Generation  of  the  megPAC  model 

Starting  from  each  cortical  MEG  source  time  series  (top  row):  1)  the  troughs  and  peaks  of 
the  optimal  low-frequency  phase  were  determined  from  the  original  sources  signal  (second 
row);  2)  the  amplitude  of  the  same  source  signal  in  the  [80,150]  Hz  gamma  range  was 
extracted  (third  row).  For  the  megPAC  model,  the  gamma  amplitude  values  at  the  troughs 
and  peaks  of  the  low-frequency  were  linearly  interpolated  (fourth  row). 
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Figure  2.  Ubiquitous  Phase-Amplitude  Coupling  in  the  Human  Cortex 

(a):  Average  low  frequency  component  for  PAC  phase,  following  anatomical  co-registration 
of  individual  cortical  surfaces  (n  =  12)  to  a  template  brain,  (b):  Representative  spectrograms 
of  ongoing  cortical  activity  averaged  at  the  phase  troughs  of  the  low  frequency  with  the 
highest  PAC  at  each  tested  brain  location  in  1  representative  subject.  The  scales  are  [-0.5, 
0.5]  s  (x-axis),  and  [40,  250]  Hz  (y-axis).  The  tested  cortical  locations  are  shown  with 
colored  dots  on  the  individual's  brain  surface,  with  color  correspondence  with  the  time- 
frequency  maps  (shown  at  lower  left  corner).  The  spectrograms  and  ongoing  cortical  signals 
were  averaged  about  the  troughs  of  the  optimal  low-frequency  oscillations  for  phase.  At 
each  cortical  location,  the  resulting  event-related  average  of  the  raw  source  signal  is  shown, 
with  the  value  of  the  corresponding  frequency  for  phase  in  the  lower  sub-panel  of  each  plot. 
The  color  bar  immediately  to  the  left  of  the  time-frequency  maps  indicates  significant 
correlation  ( p  <  0.05,  negative  in  red,  positive  in  black)  between  the  amplitude  modulations 
of  high-frequency  oscillations  and  the  event-related  average  of  the  raw  signal. 
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Figure  3.  Group  Modes  and  Patterns  of  Resting-State  Connectivity 

The  anatomy  of  the  resting-state  networks  revealed  by  the  megPAC  signal  model  based  on 
phase-amplitude  coupling  mechanisms  between  neural  ensembles  is  similar  to  those 
typically  found  using  hemodynamic  imaging  techniques.  Five  principal  spatial  modes  of 
connectivity  from  the  megPAC  signal  model  are  displayed.  For  comparison  the 
default-mode  network  (DMN),  dorsal  attention  network  combined  with  sensorimotor 
network  (DAN/SMN),  visual  network,  and  right  fronto-parietal  network  are  shown 
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based  on  fMRI  results,  as  provided  in  Vincent  et  al.,  2008  and  Yeo  et  al.,  2011.  For  the 
first  four  MEG  networks  there  is  an  overlap  to  these  four  fMRI  networks. 
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Figure  4. 

(a)  The  frontal  lobe  connections  identified  from  the  megPAC  time  series  are  consistent  with 
results  from  MRI  tractography  and  post-mortem  dissection  techniques  (Catani  et  al.,  2012). 

(b)  The  ratio  of  the  average  pairwise  signal  correlation  between  regions  within  the  DMN 
over  the  pairwise  correlation  between  brain  regions  outside  the  DMN  was  computed  for  the 
megPAC,  the  megPAC  model  with  the  low-frequency  fixed  at  5.4  Hz,  the  megPAC  model 
with  random  events  for  the  low-frequency  and  the  megPAC  model  from  random  data.  The 
regions  of  interest  within  the  DMN  are  shown  in  figure  3  and  were  selected  based  on 
Vincent  et  al.,  2008.  Overall  the  megPAC  model  yielded  the  highest  correspondence  with 
the  DMN  (*:  p  <  0.001). 
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Figure  5. 

Principles  of  Synchronized  Gating  (SG)  for  long-range  neural  communication.  The  colored 
dots  represent  neuronal  populations;  the  connecting  arrows  indicate  when  two  neuronal 
groups  can  communicate.  When  crossed  out,  this  indicates  that  neural  communication 
cannot  occur.  In  the  present  example,  only  the  red  and  blue  neuronal  groups  communicate 
with  each  other:  The  respective  phase  of  the  local  low-frequency  oscillatory  activity  are 
synchronized  and  therefore  provide  the  necessary  gating  mechanism  for  long-range 
communication  between  the  two  groups.  In  addition,  the  amplitudes  of  the  high-frequency 
bursts  are  coherent,  which  indicates  similar  levels  of  local  post-synaptic  integration  and 
spiking.  Under  the  SG  assumption,  the  green  and  black  neuronal  groups  do  not  communicate 
with  the  other  populations,  because  their  low-frequency  phase  is  not  aligned  or  the  local 
high-frequencies  are  not  coherent  with  the  red  and  blue  populations,  respectively.  In  contrast 
according  to  the  CTC  hypotheses,  the  black  and  blue  neuronal  group  would  communicate 
with  each  other. 
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